
# generate appendix C5, table C7 and figures C9 to C12

library(rio) # to load data
library(tidyverse) # to reshape data
library(lmtest) # to cluster standard errors
library(sandwich) # to cluster standard errors
library(texreg) # to generate tables
library(ggeffects) # to plot results
library(scales) # to make nice axis breaks

# set your working directory
setwd("~/replication_files")

monthly_data <- import("data/monthly_data.csv") %>%
  mutate(discovery_lag = as.factor(discovery_lag)) # to make sure that the interaction plots treat this variable as a factor

# table C7

# model 1
probit1 <- glm(sovereign_issued_1yr ~ 
               resource_rents_perc_lag*minister_is_econ_technocrat_usa + 
               log_oil_production_lag*minister_is_econ_technocrat_usa + 
               terms_of_trade_lag*minister_is_econ_technocrat_usa + 
               discovery_lag*minister_is_econ_technocrat_usa + 
               minister_is_econ_technocrat_usa +
               minister_turnover_5yrs +
               debt_crisis_experience +
               election_month + left_exec +
               fiscal_council + 
               polcon + 
               imf_agreement +
               fiscal_balance_lag + 
               tax_perc_gdp_cepal_lag +
               log_core_inflation_lag +
               gdp_pc_cepal_lag +
               gdp_growth_cepal_lag +
               ka_open_lag +
               log_reserves_lag +
               treasury_rate_lag +
               I(iso3c) + 
               time, 
               data = monthly_data,
               family = binomial(link = "probit"))

probit1_r <- probit1 %>%
  coeftest(vcovHC(probit1, type = 'HC0', cluster = c('iso3c','Summary')))

# model 2
tobit1 <- tobit(log_sovereign_issued_1yr_amount ~ 
                resource_rents_perc_lag*minister_is_econ_technocrat_usa + 
                log_oil_production_lag*minister_is_econ_technocrat_usa + 
                terms_of_trade_lag*minister_is_econ_technocrat_usa +
                discovery_lag*minister_is_econ_technocrat_usa + 
                minister_is_econ_technocrat_usa +
                minister_turnover_5yrs +
                debt_crisis_experience +
                election_month + left_exec +
                fiscal_council + 
                polcon + 
                imf_agreement +
                fiscal_balance_lag + 
                tax_perc_gdp_cepal_lag +
                log_core_inflation_lag +
                gdp_pc_cepal_lag +
                gdp_growth_cepal_lag +
                ka_open_lag +
                log_reserves_lag +
                treasury_rate_lag +
                I(iso3c) + 
                time, 
                data = monthly_data, robust = T)

# generate table C7
texreg(l = list(probit1_r,tobit1), stars = c(0.01, 0.05, 0.1),
       custom.coef.map = list("resource_rents_perc_lag" = "Resource Rents, % of GDP $_{t-1}$",
                              "log_oil_production_lag" = "Ln Oil and Gas Production $_{t-1}$",
                              "terms_of_trade_lag" = "Commodity Price Index $_{t-1}$",
                              "discovery_lag1" = "Field Discovery $_{t-1}$",
                              "minister_is_econ_technocrat_usa" = "Mainstream Minister = 1",
                              "resource_rents_perc_lag:minister_is_econ_technocrat_usa" = "Mainstream Minister X Resource Rents",
                              "minister_is_econ_technocrat_usa:log_oil_production_lag" = " Mainstream Minister X Ln Oil and Gas Production",
                              "minister_is_econ_technocrat_usa:terms_of_trade_lag" = "Mainstream Minister X Commodity Price Index",
                              "minister_is_econ_technocrat_usa:discovery_lag1" = "Mainstream Minister X Field Discovery",
                              "minister_turnover_5yrs" = "Minister Turnover (5 Years)",
                              "debt_crisis_experience" = "Debt Crisis Experience = 1",
                              "election_month" = "Election Month = 1",
                              "left_exec" = "Left Executive = 1",
                              "cb_independence" = "Central Bank Independence",
                              "fiscal_council" = "Fiscal Council = 1",
                              "polcon" = "Political Constraints",
                              "imf_agreement" = "IMF Agreement = 1",
                              "fiscal_balance_lag" = "Fiscal Balance, % of GDP $_{t-1}$",
                              "tax_perc_gdp_cepal_lag" = "Tax Revenue, % of GDP $_{t-1}$",
                              "log_core_inflation_lag" = "Ln Core Inflation $_{t-1}$",
                              "gdp_pc_cepal_lag" = "GDP Per Capita $_{t-1}$",
                              "gdp_growth_cepal_lag" = "GDP Growth, % $_{t-1}$",
                              "ka_open_lag" = "Capital Openness $_{t-1}$",
                              "log_reserves_lag" = "Ln International Reserves $_{t-1}$",
                              "Log(scale)" = "Log(Scale)"),
       fontsize = "footnotesize", digits = 3)

# retrieve fit stats
texreg(l = list(probit1,tobit1), stars = c(0.01, 0.05, 0.1))


# make predictions with clustered standard errors
# based on model 1 of table C7

# figure C9
predictions1 <- ggpredict(probit1, terms = c("resource_rents_perc_lag","minister_is_econ_technocrat_usa[all]","iso3c[all]"), vcov_fun = "vcovHC",
                          vcov_type = "HC0", 
                          vcov_args = list(cluster = monthly_data$iso3c))
plot(predictions1, grid = T) +
  labs(y = "Predicted Probability of Debt Issuance", x = bquote(`Resource Rents, % of GDP `[t-1]),
       title = " ", color = "Mainstream Minister") + theme_classic() +
  theme(axis.text = element_text(size = 12), plot.title = element_text(face = "bold"),
        axis.title = element_text(size = 12), legend.position="bottom") + 
  scale_x_continuous(breaks = pretty_breaks()) + scale_y_continuous(limits = c(0,1), labels = percent)

ggsave("figures/figureC9.pdf", height = 10, width = 8)

# figure C10
predictions2 <- ggpredict(probit1, terms = c("log_oil_production_lag","minister_is_econ_technocrat_usa[all]","iso3c[all]"), vcov_fun = "vcovHC",
                          vcov_type = "HC0", 
                          vcov_args = list(cluster = monthly_data$iso3c))
plot(predictions2, grid = T) +
  labs(y = "Predicted Probability of Debt Issuance", x = bquote(`Ln Oil and Gas Production `[t-1]),
       title = " ", color = "Mainstream Minister") + theme_classic() +
  theme(axis.text = element_text(size = 12), plot.title = element_text(face = "bold"),
        axis.title = element_text(size = 12), legend.position="bottom") + 
  scale_x_continuous(breaks = pretty_breaks()) + scale_y_continuous(limits = c(0,1), labels = percent)

ggsave("figures/figureC10.pdf", height = 10, width = 8)

# figure 11
predictions3 <- ggpredict(probit1, terms = c("terms_of_trade_lag","minister_is_econ_technocrat_usa[all]","iso3c[all]"), vcov_fun = "vcovHC",
                          vcov_type = "HC0", 
                          vcov_args = list(cluster = monthly_data$iso3c))
plot(predictions3, grid = T) +
  labs(y = "Predicted Probability of Debt Issuance", x = bquote(`Commodity Price Index `[t-1]),
       title = " ", color = "Mainstream Minister") + theme_classic() +
  theme(axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 8),
        plot.title = element_text(face = "bold"),
        axis.title = element_text(size = 12), legend.position="bottom") + 
  scale_x_continuous(breaks = pretty_breaks()) + scale_y_continuous(limits = c(0,1), labels = percent)

ggsave("figures/figureC11.pdf", height = 10, width = 8)

# figure 12
predictions4 <- ggpredict(probit1, terms = c("discovery_lag","minister_is_econ_technocrat_usa[all]","iso3c[all]"), vcov_fun = "vcovHC",
                          vcov_type = "HC0", 
                          vcov_args = list(cluster = monthly_data$iso3c))
plot(predictions4, grid = T) +
  labs(y = "Predicted Probability of Debt Issuance", x = bquote(`Field Discovery `[t-1]),
       title = " ", color = "Mainstream Minister") + theme_classic() +
  theme(axis.text = element_text(size = 12), plot.title = element_text(face = "bold"),
        axis.title = element_text(size = 12), legend.position="bottom") + 
  scale_y_continuous(limits = c(0,1), labels = percent)

ggsave("figures/figureC12.pdf", height = 10, width = 8)